﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using SmartMathLibrary.Statistics.RandomNumberGenerator;

namespace SmartMathLibrary.NumericalIntegration
{
    /// <summary>
    /// This class provides the numerical integration of multivariate integrals by using 
    /// a Monte Carlo method.
    /// </summary>
    [Serializable]
    public class MultivariateMCIntegrator
    {
        /// <summary>
        /// Internal constant.
        /// </summary>
        private const int MINPOINTS = 60;

        /// <summary>
        /// Internal constant.
        /// </summary>
        private const int MINSPLIT = 15;

        /// <summary>
        /// Internal constant.
        /// </summary>
        private const int EXPLORE_MIN = 10;

        /// <summary>
        /// Internal constant.
        /// </summary>
        private const int EXPLORE_MAX = 100;

        /// <summary>
        /// Internal constant.
        /// </summary>
        private const double EXPLORE_PROPORTION = 0.10;

        /// <summary>
        /// Internal constant.
        /// </summary>
        private const double VERY_LARGE = 1e+20;

        /// <summary>
        /// This field holds the integral.
        /// </summary>
        private MultivariateIntegral integral;

        /// <summary>
        /// This field holds the random engine.
        /// </summary>
        private readonly XORShiftRandomGenerator random;

        /// <summary>
        /// Initializes a new instance of the <see cref="MultivariateMCIntegrator"/> class.
        /// </summary>
        /// <param name="integral">The integral to integrate.</param>
        public MultivariateMCIntegrator(MultivariateIntegral integral)
        {
            if (integral == (MultivariateIntegral) null)
            {
                throw new ArgumentNullException("integral");
            }

            if (integral.A == (GeneralVector) null)
            {
                throw new ArgumentNullException("integral");
            }

            if (integral.B == (GeneralVector) null)
            {
                throw new ArgumentNullException("integral");
            }

            this.integral = integral;
            this.random = new XORShiftRandomGenerator();
        }

        /// <summary>
        /// Gets or sets the integral to integrate.
        /// </summary>
        /// <value>The integral to integrate.</value>
        public MultivariateIntegral Integral
        {
            get { return this.integral; }
            set { this.integral = value; }
        }

        /// <summary>
        /// Integrates the specified integral in a numerical way.
        /// </summary>
        /// <param name="numberOfPoints">The maximum number of points to use for the integration.</param>
        /// <returns>The value of the integral.</returns>
        public double Integrate(int numberOfPoints)
        {
            return this.RecursiveIntegration(this.integral.A, this.integral.B, this.integral.A.Count, numberOfPoints);
        }

        private double RecursiveIntegration(GeneralVector lo, GeneralVector hi, int dim, int N)
        {
            int randD = 0;
            double totalvar = 0;
            double var0 = 0;
            GeneralVector midpoint = new GeneralVector(dim);
            double[,] min = new double[2,dim];
            double[,] max = new double[2,dim];

            if (N < MINPOINTS)
            {
                return this.Integrate(lo, hi, dim, N);
            }

            SetupIntegration(lo, hi, midpoint, min, max, dim, ref N, ref randD);
            this.ExploreIntegral(midpoint, min, max, dim, randD);

            int split = ChooseSplit(min, max, dim, ref totalvar, ref var0);
            int points0 = (int) (N*var0/totalvar);

            if (points0 < MINSPLIT)
            {
                points0 = MINSPLIT;
            }

            double save = hi[split];

            hi[split] = midpoint[split];

            double result = this.RecursiveIntegration(lo, hi, dim, points0);

            hi[split] = save;
            save = lo[split];
            lo[split] = midpoint[split];
            result += RecursiveIntegration(lo, hi, dim, N - points0);
            lo[split] = save;

            return result;
        }

        private static void SetupIntegration(IGeneralVector lo, GeneralVector hi, GeneralVector midpoint, double[,] min,
                                             double[,] max, int dim, ref int n, ref int randD)
        {
            for (int j = 0; j < dim; j++)
            {
                midpoint[j] = 0.5*(lo[j] + hi[j]);
                min[0, j] = min[1, j] = VERY_LARGE;
                max[0, j] = max[1, j] = -VERY_LARGE;
            }

            randD = (int) (n*EXPLORE_PROPORTION);

            if (randD < EXPLORE_MIN)
            {
                randD = EXPLORE_MIN;
            }
            if (randD > EXPLORE_MAX)
            {
                randD = EXPLORE_MAX;
            }

            n -= randD;
        }

        private void ExploreIntegral(GeneralVector midpoint, double[,] min, double[,] max, int dim, int numberOfSamples)
        {
            GeneralVector point = new GeneralVector(dim);

            for (int n = 0; n < numberOfSamples; n++)
            {
                for (int j = 0; j < dim; j++)
                {
                    point[j] = this.integral.A[j] + this.random.NextDouble()*(this.integral.B[j] - this.integral.A[j]);
                }

                double fx = this.integral.Function.SolveAt(point);

                for (int j = 0; j < dim; j++)
                {
                    int half = 0;

                    if (point[j] > midpoint[j])
                    {
                        half = 1;
                    }

                    max[half, j] = fx > max[half, j] ? fx : max[half, j];
                    min[half, j] = fx < min[half, j] ? fx : min[half, j];
                }
            }
        }

        private int ChooseSplit(double[,] min, double[,] max, int dim, ref double var, ref double var0)
        {
            int split = -1;

            for (int j = 0; j < dim; j++)
            {
                if (max[0, j] > min[0, j] && max[1, j] > min[1, j])
                {
                    double sigma0 = Math.Pow(max[0, j] - min[0, j], 2.0/3.0);
                    double sigma1 = Math.Pow(max[1, j] - min[1, j], 2.0/3.0);
                    double sigma = sigma0 + sigma1;

                    if ((split == -1) || (sigma < var))
                    {
                        split = j;
                        var = sigma;
                        var0 = sigma0;
                    }
                }
            }

            if (split == -1)
            {
                var0 = 1.0;
                var = 2.0;
                split = (int) (this.random.NextDouble()*dim);
            }

            return split;
        }

        private double Integrate(GeneralVector lo, GeneralVector hi, int dim, double n)
        {
            double sum = 0.0;
            GeneralVector point = new GeneralVector(dim);

            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < dim; j++)
                {
                    point[j] = this.integral.A[j] + random.NextDouble()*(this.integral.B[j] - this.integral.A[j]);
                }

                double fx = this.integral.Function.SolveAt(point);

                sum += fx;
            }

            double volume = 1.0;

            for (int i = 0; i < dim; i++)
            {
                volume *= (hi[i] - lo[i]);
            }

            return volume*sum/n;
        }
    }
}